clear; close all; clc
warning off;
addpath('../LibForSerial/LibrarySLP','../LibForSerial/LibraryLeSage/','../LibForSerial/LibraryRamey')

%----------- Load shocks using two identification schemes ----------------%
   load('../Figure2_4_JLN/DataShocks/dataJMNupdated.mat');  
 uF2015        = ehat_maxG(:,3); clearvars -except uF2015 datem         projectdir projectdirDB
 
%----------- Construct indicators for seq. of shocks with 2 ident schemes-%
IndPreviousNShocks       = zeros(size(uF2015));
  NumberConsecutiveshocks = 2; % 1st previous shocks that is  positive, and the 2nd one is positive too --> WP2015 P=2,6 with and without lagmatrix( data.*repmat(IndSeqShocks,1,size(data,2)) , 1:P )
%/////////////////////////////////////////////////////////////////////////
for tt=NumberConsecutiveshocks:length(uF2015)-1 
    if tt==NumberConsecutiveshocks
        if (                                        sum(uF2015(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                  );  IndPreviousNShocks(tt) = 1; end
    else
        if (uF2015(tt-NumberConsecutiveshocks)<0 && sum(uF2015(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                  );  IndPreviousNShocks(tt) = 1; end
    end
end
clear tt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sample=[ 1961 1; 2019  12];
data    = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
dataInd = xlsread('../Dataset/VARdataBBupdatedMonthly.xlsx','Indicators');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data    = data(   startsample:endsample,2:end);
dataInd = dataInd(startsample:endsample,2:end);
clear startsample endsample 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:10]) = log(data(:,[3:10]))*100;
%-------------------------------------------------------------------------%
    Remove_Large_Shocks = 0;

    IndSeqShocks = IndPreviousNShocks;        % JLN Financial Uncertainty shocks (from trivariate system with Macro Uncertainty)
    IndLargeShocks = uF2015>2;

LargeShocksToBeRemoved = IndSeqShocks & IndLargeShocks; 
if Remove_Large_Shocks;  IndSeqShocks = IndSeqShocks - LargeShocksToBeRemoved; end
%-------------------------------------------------------------------------%
T = size(data,1);

yyyy=repmat((1961:2019)',1,12);
yyyy=reshape(yyyy' ,size(yyyy,1)*12,1);
mm  =repmat((1:12)',size(yyyy,1)/12,1);
datem=yyyy*100+mm;
datem=datem(1:end);





clear yyyy mm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT SPECIFICATION CHOICES 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
timesample =1;    % 1 means 1961-2019; 2 means 1986-2019
statechoice=1;    % 1 means NBER (or Recession Probability); 2 means PMI; 3 means unemp threshold; 4 means ZLB threshold
shockchoice=1;    % 1 means JMN shock from Xt = (UM ; ip; UF) 

 trend=0;  % 4 means quartic; 2 means quadratic; 0 means no trend  
 nlag =3;  % number of lags of control variables
 opt  =0;  % 0 means Newey-West with bandwidth=i; 1 means optimal bandwidth (takes much longer to run)'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hor   = 12*3;   % horizon for IRFs
%clevel= 1.960; % confidence level 1.96 --> 95% (CI=0.95; clevel = norminv(CI + (1-CI)/2))  
 clevel= 1.645; % confidence level 1.64 --> 90% (CI=0.90; clevel = norminv(CI + (1-CI)/2))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%IMPORTANT INPUT (sample window)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if timesample==1
    t1=find(datem==sample(1,1)*100+sample(1,2)); 
    t2=find(datem==sample(2,1)*100+sample(2,2));
end


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %DATA
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vdata=data(t1:t2,:); clear data 
  plottype = 'Output';      ii=3; PlotLabel = 'Industrial Production';% 
% plottype = 'riskfr';      ii=1; %
%                          %ii=2; PlotLabel = 'Short rate';% for Wu-Xia shadow rate  
% plottype = 'Inflat';     %ii=8; PlotLabel = 'Inflation'; % CPI
%                           ii=9; PlotLabel = 'Inflation'; % PCE
% plottype = 'SP500' ;      ii=6; PlotLabel = 'Stock Market';
%                          %ii=7; PlotLabel = 'Stock Market'; % real SP500 

 tb3     =vdata(:,1);
 ipg     =vdata(:,3);
%tb3     =vdata(:,2); % for Wu-Xia shadow rate
 inflat  =vdata(:,9);
 mkt     =vdata(:,6);
%mkt     =vdata(:,7); % real SP500 
 empl    =vdata(:,10);

recession     = dataInd(:,1);
recessionProb = dataInd(:,2);
pmi           = dataInd(:,3);
unemp         = dataInd(:,4);
unemphp       = dataInd(:,5);
zlbdummy      = dataInd(:,6);% the ZLB episode (2009Q1-2015Q4) as in Debortoli, Gali, Gambetti NBER Macroeconomics Annual 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Creating variables of interest
inflat  = [NaN(12,1); inflat(13:end)-inflat(1:end-12)]; 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %IMPORTANT INPUT (threshold for bad state indicator, e.g., the unemployment rate)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   statechoice= 1;
if statechoice==1
    state=recession;  threshold=1;   %NBER (or Recession Probability) For Figure 3
elseif statechoice==2
    state=44.5;       threshold=pmi;   
elseif statechoice==3
%   state=unemp;      threshold=6.5;
    state=unemp;      threshold=unemphp;   
end

rrr=find(state>=threshold); 
%rrrn=rrr  ; %Since STATE above threshold in the current  period is the criterion
 rrrn=rrr+1; %Since STATE above threshold in the previous period is the criterion
if rrrn(end)==length(state)+1
    rrrn=rrrn(1:end-1);
else
    rrrn=rrrn;
end
fu=zeros(size(ipg));
fu(rrrn)=1; % This is the dummy variable that defines the two states

disp('Percentage of the sample assigned to Recessions')
disp((length(rrr)/length(state))*100)
clear state

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ii==1
    X1=[ipg tb3];
elseif ii==3
   %X1=[ipg tb3];
    X1=[ipg tb3 empl];
elseif ii ==6 || ii ==7
    X1=[mkt tb3];
end

[rx,cx]=size(X1); 




lags=0:1:nlag;
XLAG1 = lagmatrix(X1,[lags]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SHOCK IDENTIFICATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if shockchoice==1      %JMN AEJ shock
    shock         = uF2015(nlag+1:end);    
    disp('LMN AEJ shock (from trivariate system with IP )')
end
IndSeqShocks = IndSeqShocks(nlag+1:end);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PRELIMINARIES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

constant= ones(length(shock),1);
t=1:1:length(constant);
tsq=t.^2;
tcu=t.^3;
tqu=t.^4;

xorig= XLAG1(nlag+1:end, cx+1:end);%only picks lagged values



%Y variable
if ii==1
    yy=[tb3(nlag+1:end)]; 
elseif ii==3
    yy=[ipg(nlag+1:end)]; 
elseif ii==6 || ii ==7
    yy=[mkt(nlag+1:end)];
end

smoothLP = 1; % 0 for classical LP; 1 for BB ReStat (for paper)
r=2;
lambda=5;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LINEAR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if trend==4
    x=[constant, t', tsq', tcu', tqu', shock,  xorig]; %complete RHS
    rpos=6; % position of shock
elseif trend==2
    x=[constant, t', tsq',             shock,  xorig];
    rpos=4;
elseif trend==0
    x=[constant,                       shock,  xorig];
    rpos=2;
end
[liny, confidencey]=linlp_rz(yy,x,hor,rpos,1, clevel, opt);



fu=fu(nlag+1:end);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NON-LINEAR
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x;
x=xorig;
if trend==4 % allow quartic trend
    x=[t', tsq', tcu', tqu', constant, (1-fu).*constant, fu.*shock, (1-fu).*shock, repmat(fu,1,size(x,2)).*x, repmat((1-fu),1,size(x,2)).*x];
    rpos=7; 
elseif trend==2 % allow quadratic trend
    x=[t', tsq',             constant, (1-fu).*constant, fu.*shock, (1-fu).*shock, repmat(fu,1,size(x,2)).*x, repmat((1-fu),1,size(x,2)).*x];
    rpos=5; 
elseif trend==0 %no trend
    x=[                      constant, (1-fu).*constant, fu.*shock, (1-fu).*shock, repmat(fu,1,size(x,2)).*x, repmat((1-fu),1,size(x,2)).*x];
    rpos=3; % position of shock
end
[stateay, stateby, confidenceya, confidenceyb]=statelp_rz(yy,x,hor,rpos,1, clevel, opt); 
if smoothLP
    slp_b    = locproj(        yy,x(:,rpos),            x(:,setdiff(2:size(x,2),[rpos]))  ,0,hor,'smooth',r,lambda); %IR from Smooth Local Projection 
    slp_g    = locproj(        yy,x(:,rpos+1),          x(:,setdiff(2:size(x,2),[rpos+1])),0,hor,'smooth',r,lambda); %IR from Smooth Local Projection 
%plot([stateay'   slp_b.IR(1:end-1)])
       stateay = (slp_b.IR(1:end-1))';
%plot([confidenceya'  [slp_b.Interval_down(1:end-1), slp_b.Interval_up(1:end-1)]])       
       confidenceya = [slp_b.Interval_down(1:end-1), slp_b.Interval_up(1:end-1)]';
%plot([stateby'   slp_g.IR(1:end-1)])
       stateby = (slp_g.IR(1:end-1))';
%plot([confidenceyb'  [slp_g.Interval_down(1:end-1), slp_g.Interval_up(1:end-1)]])       
       confidenceyb = [slp_g.Interval_down(1:end-1), slp_g.Interval_up(1:end-1)]';    
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% NON-LINEAR and Consecutive Shocks
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear x;
x=xorig;
x=[constant, (1-fu).*constant, fu.*shock, fu.*shock.*IndSeqShocks, (1-fu).*shock, (1-fu).*shock.*IndSeqShocks, repmat(fu,1,size(x,2)).*x, repmat((1-fu),1,size(x,2)).*x];
rpos=3; % position of shock

[stateay_cons, stateby_cons, confidenceya_cons, confidenceyb_cons, stateayMult, confidenceyaMult, stateay_bis]=statelp_rz_consec(yy,x,hor,rpos,1, clevel, opt); 
if smoothLP
slp_sd_b = locprojStateDep(yy,x(:,rpos)  ,x(:,rpos+1),x(:,setdiff(2:size(x,2),[rpos   rpos+1])),0,hor,'smooth',r,lambda); %IR from Smooth Local Projection
slp_sd_g = locprojStateDep(yy,x(:,rpos+2),x(:,rpos+3),x(:,setdiff(2:size(x,2),[rpos+2 rpos+3])),0,hor,'smooth',r,lambda); %IR from Smooth Local Projection
%plot([stateay_cons'        slp_sd_b.IR(1:end-1)+slp_sd_b.IRstate(1:end-1)])
       stateay_cons      = (slp_sd_b.IR(1:end-1)+slp_sd_b.IRstate(1:end-1))';
       stateayMult       =                       slp_sd_b.IRstate(1:end-1)' ;
%plot([confidenceya_cons'  [slp_sd_b.Interval_down(1:end-1)+slp_sd_b.Interval_down_state(1:end-1), slp_sd_b.Interval_up(1:end-1)+slp_sd_b.Interval_up_state(1:end-1)]])       
       confidenceya_cons = [slp_sd_b.Interval_down(1:end-1)+slp_sd_b.Interval_down_state(1:end-1), slp_sd_b.Interval_up(1:end-1)+slp_sd_b.Interval_up_state(1:end-1)]';
       confidenceyaMult  = [                                slp_sd_b.Interval_down_state(1:end-1),                               slp_sd_b.Interval_up_state(1:end-1)]';
%plot([stateby_cons'        slp_sd_g.IR(1:end-1)+slp_sd_g.IRstate(1:end-1)])
       stateby_cons      = (slp_sd_g.IR(1:end-1)+slp_sd_g.IRstate(1:end-1))';
%plot([confidenceyb_cons'  [slp_sd_g.Interval_down(1:end-1)+slp_sd_g.Interval_down_state(1:end-1), slp_sd_g.Interval_up(1:end-1)+slp_sd_g.Interval_up_state(1:end-1)]])       
       confidenceyb_cons = [slp_sd_g.Interval_down(1:end-1)+slp_sd_g.Interval_down_state(1:end-1), slp_sd_g.Interval_up(1:end-1)+slp_sd_g.Interval_up_state(1:end-1)]';
end       

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zz=zeros(1,hor);
i=1;

%% Figure 3, Panels (a), (b), and (c)
FontSize = 46;
LineWidth =  8;
MarkerSize= 30;

figure(31); hold on
 h1=line(1:1:hor, stateay(i,:), 'color', [0 0 0]+0.05*12,'LineStyle','-','linewidth',LineWidth);
 h2=line(1:1:hor, liny(i,:), 'color','k'        ,'LineStyle','-' ,'linewidth',LineWidth);
 h3=line(1:1:hor, stateby(i,:), 'color', [0 0 0]+0.05*12,'LineStyle',':','linewidth',LineWidth);
hold on; plot(1:hor , zz , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
ylabel(PlotLabel,'FontSize',FontSize);
xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
                   set(gca,                  'FontSize',FontSize)
% For paper
  set(gca,'FontSize',52)
  hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(1),'MarkerSize',32); set(hline(3),'MarkerSize',32)
% hAxes = gca;
% reformatTickLabels(hAxes);
%[legh,objh] = legend([h2 h1 h3],{'Linear','Bad times','Good times'},'Fontsize',30,'Location','Best');
 [legh,objh] = legend([h2 h1 h3], 'Linear','Bad times','Good times',               'Location','Best');
 set(legh,'FontSize',36);
lineh = findobj(objh,'type','line');
set(lineh,'LineWidth',3)
%ylim([-2 0.4]) for plot D.1

figure(32)
 line(1:1:hor, stateay(     i,:), 'color', [0 0 0]+0.05*12,'LineStyle' ,'-','linewidth',LineWidth)
 line(1:1:hor, stateay_cons(i,:), 'color', [0 0 0]+0.05*12,'LineStyle','--','linewidth',LineWidth)
hold on; plot(1:hor , zz , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
%ylabel(PlotLabel,'FontSize',FontSize);
xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
                   set(gca,                  'FontSize',FontSize)
% For paper
  set(gca,'FontSize',52)
  hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(1),'MarkerSize',32); set(hline(2),'MarkerSize',32)
% hAxes = gca;
% reformatTickLabels(hAxes);
%[legh,objh] = legend({'Bad times','Bad times (two consec. shocks)'},'Fontsize',30,'Location','Best');
 [legh,objh] = legend( 'Bad times','Bad times (two consec. shocks)',               'Location','Best');
 set(legh,'FontSize',36);
lineh = findobj(objh,'type','line');
set(lineh,'LineWidth',3)

figure(33)
grpyat=[(1:1:hor)', confidenceyaMult(1,:,i)'; (hor:-1:1)' confidenceyaMult(2,hor:-1:1,i)']; patch(grpyat(:,1), grpyat(:,2), [0.7 0.7 0.7],'edgecolor', [0.7 0.7 0.7]);
hold on
h1= line(1:1:hor, stateayMult(i,:),'LineStyle','--','linewidth',LineWidth,'color',[0.72 0.27 1.0]);
hold on; plot(1:hor , zz , 'k-' , 'LineWidth' , 2 ); grid; xlim([1 hor]); 
xlim([1      hor]);set(gca,'XTick',[6:6:hor],'FontSize',FontSize);
                   set(gca,                  'FontSize',FontSize)
% For paper
  set(gca,'FontSize',52)
  hline = findobj(gcf, 'type', 'line'); set(hline,'LineWidth',6); set(hline(2),'MarkerSize',32); 